Ejemplos de regresión gamma
Oihane Álvarez, Gabriel Carbonell, Daniel Hernández, Celia Sifre
08/05/2022
Análisis de la pelea final de Los Vengadores: Endgame.
Introducción
Para la realización de este ejercicio he empleado la base de datos avengers.csv de la librería flexplot. He modificado la base de datos para solo tener en cuenta las variables que se emplean en este video y les he cambiado el nombre. Si queréis la base de datos original, la podeís encontrar en este enlace.
Esta base de datos contiene originalmente los atributos de combate de 812 luchadores en la pelea final de Los Vengadores: Endgame. Para este ejercicio solo se emplearán 3 atributos, que son los siguientes: minutes.fighting (renombrado a T_Batalla), son los minutos que un luchador es capaz de durar en la batalla hasta que muere o se rinde; willpower (renombrado a Determinacion), que según la documentación, está determinado por el periódo de tiempo que el luchador es capaz de esperar en la DMV (aquí en España es la DGT) para que le entreguen el carnet de conducir; y finalmente, injuries (renombrado a Heridas) que es el número de heridas que pueden soportar.
Escena de la batalla final de Los Vengadores: Endgame
En este ejercicio vamos a tratar de modelizar el tiempo en batalla (T_batalla) en función de las variables Determinación y Heridas empleando la regresión gamma.
En primer lugar vamos a cargar el banco de datos, y siguiendo la forma en que trata a la variable Heridas en el vídeo del enlace del primer párrafo, vamos a categorizar en tres grupos la variable Heridas mediante la función cut().
datos <- read.csv("avengers_mod.csv");datos <- datos[-437,] ## valor negativo de determinación
datos$Heridas <- cut(datos$Heridas, breaks = c(-1,2,4,5),labels = c("0-2","3-4","5+"))Con esto realizado, podemos empezar con la parte de estadística descriptiva del banco de datos.
Análisis previo descriptivo del banco de datos
En primer lugar vamos a obtener un resumen numérico de las variables mediante la función summary().
summary(datos)## T_batalla Determinacion Heridas
## Min. : 1.10 Min. : 4.00 0-2:206
## 1st Qu.: 8.70 1st Qu.: 48.00 3-4:374
## Median :10.40 Median : 60.00 5+ :231
## Mean :12.38 Mean : 60.07
## 3rd Qu.:13.15 3rd Qu.: 72.00
## Max. :93.70 Max. :115.00Observamos que no hay ningún valor ausente en ninguna de las variables. La mayoría de los valores de nuestra variable respuesta T_batalla se concentran en el intervalo [1.10,13.12] quedando algunos de estos valores muy alejados del resto. En cuanto a los valores de la variable Determianción podemos intuir que son bastante simétricos, con media en 59.99. Finalmente, en cuanto a la variable Heridas podemos decir que la hay un número similar de luchadores con 0-2 heridas, que con 3-4 y que con 5+,siendo el grupo con 3-4 heridas el más numeroso.
Pasamos ahora a la parte gráfica del análisis descriptivo. Vamos a comenzar esta parte obteniendo el histograma T_batalla.
ggplot(data = datos, aes(T_batalla)) +
geom_histogram(color = "dodgerblue", fill = "skyblue") +
labs(title = "Figura 1. Histograma de T_batalla") +
theme_bw()
El histograma cumple con una de las particularidades de la distribución gamma: es asimétrica por la derecha. Observamos que las probabilidades de que un luchador esté mucho tiempo peleando son muy bajas, durando la mayoría de ellos tiempos comprendidos entre los 0 y los 25 minutos.
Ahora vamos a ver qué tipo de relación hay entre las variables. Para ello, vamos a obtener el diagrama de dispersión entre T_batalla y Determinacion .
ggplot(data = datos, aes(x = Determinacion, y = T_batalla)) +
labs(title = "Figura 2. Diagrama de dispersión T_batalla frente a Determinacion") +
geom_point() + theme_bw()
A través de este diagrama podemos intuir que cuanto mayor es la Determinacion mayor es el tiempo en batalla. Sin embargo, no está del todo claro. Así que vamos a ver si hay alguna relación si añadimos la variable Heridas y obtenemos un diagrama de dispersión para cada grupo.
ggplot(data = datos, aes(x = Determinacion, y = T_batalla)) +
labs(title = "Figura 3. Diagramas de dispersión de T_batalla frente a Determinacion por grupos de Heridas") +
geom_point() + facet_wrap(~Heridas) + theme_bw() +
theme(plot.title = element_text(size=11))
Ahora la relación parece bastante más clara. A mayor Determinación máyor es T_batalla, solo que la variable Heridas disminuye el efecto que tiene la Determinacion, de modo que a mayor número de heridas, menos importante es el efecto de Determinacion. Es decir, parece que hay interacción entre las variables Determinacion y Heridas.
Modelización Lineal
Aunque sepamos que estamos analizando una variable continua no negativa, ajustamos los datos a un modelo lineal para observar si se ajusta adecuadamente.
mod.lineal <- lm(T_batalla ~ ., data=datos)
summary(mod.lineal)##
## Call:
## lm(formula = T_batalla ~ ., data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.808 -2.913 -0.711 1.286 75.261
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.49952 1.07994 10.648 < 2e-16 ***
## Determinacion 0.10676 0.01469 7.269 8.55e-13 ***
## Heridas3-4 -6.33453 0.64696 -9.791 < 2e-16 ***
## Heridas5+ -9.15367 0.71184 -12.859 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.376 on 807 degrees of freedom
## Multiple R-squared: 0.2411, Adjusted R-squared: 0.2383
## F-statistic: 85.48 on 3 and 807 DF, p-value: < 2.2e-16par(mfrow=c(2,2)); plot(mod.lineal);par(mfrow=c(1,1))
Como podemos observar en los plots de diagnóstico, los residuos no se ajustan adecuadamente. En primer lugar, el gráfico residuals vs. fitted no sugiere una relación lineal perfecta, puesto que hay residuos que no siguen la línea horizontal. muchos de los residuos se distribuyen a lo largo de la línea horizontal. Por otra parte, el QQ plot muestra como la mayoría de los residuos caen en la línea de referencia, pero no así el extremo superior. Esto sugiere que los datos siguen una distribución de cola larga. Tampoco parece que se cumpla adecuadamente la hipótesis de homocedasticidad. El gráfico del leverage muestra también muchos valores extremos. Por tanto, la regresión gaussiana no parece la mejor alternativa para modelizar estos datos.
Modelización con Regresión Gamma
Como estamos analizando una variable continua no negativa en la que de forma general su varianza parece aumentar con la media, parece razonable emplear un ajuste con una regresión gamma.
La variable respuesta \(Y_{ij}\) representa el número de minutos que aguanta en combate el luchador i perteneciente al grupo j de la variable Heridas.
La variable respuesta \(Y_{ij}\) se distribuye Gamma de parámetros \(\nu\) y \(\lambda\):
\[ Y_{ij} \sim Ga(\nu_{ij}, \lambda_{ij}) \]
La componente sistemática del modelo está formada por 2 variables explicativas (Determinación, variable cuantitativa continua, y Heridas, variable categórica ordinal).
Sabemos que \(\lambda = \nu / \mu\), por lo que \(E(Y_{ij})=\frac{\nu}{\lambda}=\mu_{ij}\). Para relacionar el predictor lineal con la respuesta media, podemos emplear 3 funciones link:
- La función inversa: \(g(\mu_{ij}) = \frac{1}{\mu_{ij}}\)
- La función log: \(g(\mu_{ij}) = \log \mu_{ij}\)
- La función identidad: \(g(\mu_{ij}) = \mu_{ij}\)
De tal forma que el predictor lineal queda: \[ g(\mu_{ij}) = \beta_0 + \beta_1 x_{i} + \gamma_1 d_j +\gamma_2 d_j \]
donde \(x_{ij}\) es la Determinación del individuo i y dj es una variable indicadora para la que \(\gamma_1 \cdot 1\) si Heridas = 3-4 y \(\gamma_2 \cdot1\) si Heridas = 5+.
Ajuste con función link inversa
Como el modelo no tiene muchas covariables, vamos a ajustar manualmente 4 modelos distintos: un primer modelo reducido en el que solo hay una variable explicativa (Determinacion); un segundo modelo completo con solo los efectos pricipales de las dos variables explicativas; un tercer modelo con los efectos principales y la interacción entre las variables explicativas; y un cuarto modelo con el efecto principal de la variable Determinacion y su interacción con la variable Heridas.
inv.reducido <- glm(T_batalla ~ Determinacion, data = datos, family = Gamma(link = "inverse"))
summary(inv.reducido)##
## Call:
## glm(formula = T_batalla ~ Determinacion, family = Gamma(link = "inverse"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.65008 -0.30653 -0.12254 0.06745 2.98019
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.331e-01 6.398e-03 20.812 <2e-16 ***
## Determinacion -8.254e-04 9.172e-05 -8.998 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.37993)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 157.05 on 809 degrees of freedom
## AIC: 4898
##
## Number of Fisher Scoring iterations: 6par(mfrow = c(2,2), mar = c(2.1,4,2.1,4))
plot(inv.reducido)
inv.completo <- glm(T_batalla ~ Determinacion + Heridas,
data = datos, family = Gamma(link = "inverse"))
summary(inv.completo)##
## Call:
## glm(formula = T_batalla ~ Determinacion + Heridas, family = Gamma(link = "inverse"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.49902 -0.25253 -0.06299 0.11087 2.77437
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.635e-02 5.147e-03 18.720 <2e-16 ***
## Determinacion -6.143e-04 6.708e-05 -9.157 <2e-16 ***
## Heridas3-4 2.864e-02 2.750e-03 10.414 <2e-16 ***
## Heridas5+ 5.663e-02 3.927e-03 14.421 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2102858)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 102.87 on 807 degrees of freedom
## AIC: 4549.9
##
## Number of Fisher Scoring iterations: 6inv.interac <- glm(T_batalla ~ Determinacion*Heridas,
data = datos, family = Gamma(link = "inverse"))
summary(inv.interac)##
## Call:
## glm(formula = T_batalla ~ Determinacion * Heridas, family = Gamma(link = "inverse"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.48389 -0.24876 -0.06814 0.11584 2.72977
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.849e-02 6.832e-03 14.415 < 2e-16 ***
## Determinacion -6.437e-04 9.084e-05 -7.086 3.02e-12 ***
## Heridas3-4 1.939e-02 1.007e-02 1.926 0.0545 .
## Heridas5+ 6.683e-02 1.444e-02 4.627 4.32e-06 ***
## Determinacion:Heridas3-4 1.434e-04 1.464e-04 0.980 0.3274
## Determinacion:Heridas5+ -1.668e-04 2.133e-04 -0.782 0.4345
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2087984)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 102.43 on 805 degrees of freedom
## AIC: 4550.3
##
## Number of Fisher Scoring iterations: 6inv.red.int <- glm(T_batalla ~ Determinacion + Determinacion:Heridas,
data = datos, family = Gamma(link = "inverse"))
summary(inv.red.int)##
## Call:
## glm(formula = T_batalla ~ Determinacion + Determinacion:Heridas,
## family = Gamma(link = "inverse"), data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.55654 -0.24644 -0.07537 0.10053 2.71392
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.157e-01 4.727e-03 24.484 <2e-16 ***
## Determinacion -8.608e-04 6.199e-05 -13.887 <2e-16 ***
## Determinacion:Heridas3-4 3.927e-04 4.148e-05 9.466 <2e-16 ***
## Determinacion:Heridas5+ 7.956e-04 6.175e-05 12.885 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2198918)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 107.10 on 807 degrees of freedom
## AIC: 4583.2
##
## Number of Fisher Scoring iterations: 6En general, la aplicabilidad de estos modelos falla en la hipótesis de normalidad. Vemos que el diagrama de cuantiles normales indica que los residuos grandes se desvían mucho de la línea teórica, y a simple vista parece que hay problemas leves de homocedasticidad, aunque no muy preocupantes. Por temas de espacio, solo hemos incluído las gráficas del primer modelo, las del resto de modelos son muy similares, así que no aporta mucho para el espacio que ocupan.
Ahora vamos a probar a introducir términos polinómicos de la variable Determinacion con tal de ver si mejora tanto el ajuste como la aplicabilidad. Vamos a automatizar la búsqueda del mejor modelo con términos polinómicos, de tal forma que nos quedaremos con aquel que tenga menor AIC.
AICs <- c()
for(i in 1:6){
model <- glm(T_batalla ~ poly(Determinacion,i), data = datos,
family = Gamma(link = "inverse"))
AICs[i] <- model$aic
}
AICs## [1] 4898.008 4898.219 4900.188 4901.182 4896.510 4897.431Atendiendo a los valores del AIC el mejor modelo sería el del polinomio de grado 6, sin embargo, la diferencia con el modelo de grado 1 es insignificante, así que por simplicidad, nos quedamos con el modelo del polinomio de grado 1.
Vamos a probar también con la transformación logaritmica de Determinación y ver si mejora el modelo.
summary(mod.log <- glm(T_batalla ~ log(Determinacion),
data = datos, family = Gamma(link = "inverse")))##
## Call:
## glm(formula = T_batalla ~ log(Determinacion), family = Gamma(link = "inverse"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.65494 -0.31558 -0.13255 0.07621 2.91980
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.282949 0.024278 11.65 <2e-16 ***
## log(Determinacion) -0.049265 0.005823 -8.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.3753669)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 157.35 on 809 degrees of freedom
## AIC: 4899.6
##
## Number of Fisher Scoring iterations: 6Estrictamente hablando la transformación logaritmica de Determinación mejora el ajuste del modelo, así que en este caso vamos a conservarla porque no añade complejidad al modelo.
A continuación vamos a introducir la variable Heridas al modelo con la transformación logaritmica. Vamos a ajustar un primer modelo con solo los efectos principales y un segundo modelo que incluye también la interacción.
summary(mod.log.completo <- glm(T_batalla ~ log(Determinacion)+Heridas,
data = datos, family = Gamma(link = "inverse")))##
## Call:
## glm(formula = T_batalla ~ log(Determinacion) + Heridas, family = Gamma(link = "inverse"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.50089 -0.25688 -0.07034 0.11060 2.77487
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.210178 0.017976 11.692 <2e-16 ***
## log(Determinacion) -0.037260 0.004227 -8.815 <2e-16 ***
## Heridas3-4 0.028810 0.002718 10.599 <2e-16 ***
## Heridas5+ 0.056799 0.003885 14.619 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2067435)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 102.55 on 807 degrees of freedom
## AIC: 4547.2
##
## Number of Fisher Scoring iterations: 6summary(mod.log.inter <- glm(T_batalla ~ log(Determinacion)*Heridas,
data = datos, family = Gamma(link = "inverse")))##
## Call:
## glm(formula = T_batalla ~ log(Determinacion) * Heridas, family = Gamma(link = "inverse"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.49015 -0.24910 -0.07004 0.10942 2.71231
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.235762 0.026753 8.812 < 2e-16 ***
## log(Determinacion) -0.043292 0.006283 -6.890 1.12e-11 ***
## Heridas3-4 -0.036312 0.037683 -0.964 0.3355
## Heridas5+ 0.066892 0.053165 1.258 0.2087
## log(Determinacion):Heridas3-4 0.015706 0.009034 1.739 0.0825 .
## log(Determinacion):Heridas5+ -0.002677 0.012811 -0.209 0.8346
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2037018)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 101.79 on 805 degrees of freedom
## AIC: 4545.1
##
## Number of Fisher Scoring iterations: 6De nuevo, por simplicidad, vamos a escoger el modelo sin interacciones, pese a que su AIC sea un par de unidades mayor que el que incluye las interacciones.
Finalmente, mediante validación cruzada vamos a ver qué modelo predice mejor. Dicho modelo será con el que nos quedemos. Vamos a comparar los modelos completos con y sin transformación logarítmica, pues tienen un AIC similar, si bien el de la transformación logarítmica es algo menor.
set.seed(1)
cv.glm(datos, inv.completo, K = 10)$delta## [1] 53.67619 53.63140cv.glm(datos, mod.log.completo, K = 10)$delta## [1] 53.44267 53.40633Ambos modelos cometen un error de predicción prácticamente idéntico de modo que por interpretabilidad escogemos el modelo sin la transformación logarítmica. Determinamos que proporción de la DEVIANCE explica el modelo.
1-(inv.completo$deviance/inv.completo$null.deviance)## [1] 0.446797Por lo tanto, para este modelo en el cual se ha empleado el enlace inverse, se consigue un AIC=4549 y un porcentaje de la DEVIANCE explicada del 44.67%.
Ajuste con función link log
Para este ajuste vamos a seguir un proceso idéntico al seguido con la función link inversa. Comenzamos ajustando 4 modelos: un modelo reducido, un modelo completo, un modelo con las interacciones y un modelo reducido con interacción.
log.reducido <- glm(T_batalla ~ Determinacion, data = datos, family = Gamma(link = "log"))
summary(log.reducido)##
## Call:
## glm(formula = T_batalla ~ Determinacion, family = Gamma(link = "log"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.64873 -0.31377 -0.12936 0.06931 2.94025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.862709 0.075601 24.639 <2e-16 ***
## Determinacion 0.010583 0.001206 8.772 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.3757036)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 156.73 on 809 degrees of freedom
## AIC: 4896.3
##
## Number of Fisher Scoring iterations: 4par(mfrow = c(2,2), mar = c(2.1,4,2.1,4))
plot(inv.reducido)
log.completo <- glm(T_batalla ~ Determinacion + Heridas, data = datos, family = Gamma(link = "log"))
summary(log.completo)##
## Call:
## glm(formula = T_batalla ~ Determinacion + Heridas, family = Gamma(link = "log"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.47796 -0.25894 -0.07028 0.10619 2.79336
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.374921 0.067720 35.070 <2e-16 ***
## Determinacion 0.007980 0.000921 8.664 <2e-16 ***
## Heridas3-4 -0.411992 0.040569 -10.155 <2e-16 ***
## Heridas5+ -0.696779 0.044637 -15.610 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2139531)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 103.82 on 807 degrees of freedom
## AIC: 4557.5
##
## Number of Fisher Scoring iterations: 5log.interac <- glm(T_batalla ~ Determinacion*Heridas, data = datos, family = Gamma(link = "log"))
summary(log.interac)##
## Call:
## glm(formula = T_batalla ~ Determinacion * Heridas, family = Gamma(link = "log"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.48401 -0.24846 -0.07244 0.11820 2.73054
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.060095 0.118300 17.414 < 2e-16 ***
## Determinacion 0.012792 0.001763 7.255 9.43e-13 ***
## Heridas3-4 0.029201 0.143569 0.203 0.83888
## Heridas5+ -0.334047 0.157389 -2.122 0.03411 *
## Determinacion:Heridas3-4 -0.006995 0.002213 -3.161 0.00163 **
## Determinacion:Heridas5+ -0.005627 0.002440 -2.306 0.02136 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2047699)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 101.75 on 805 degrees of freedom
## AIC: 4544.8
##
## Number of Fisher Scoring iterations: 5log.red.int <- glm(T_batalla ~ Determinacion + Determinacion:Heridas,
data = datos, family = Gamma(link = "log"))
summary(log.red.int)##
## Call:
## glm(formula = T_batalla ~ Determinacion + Determinacion:Heridas,
## family = Gamma(link = "log"), data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.52599 -0.25831 -0.07079 0.11178 2.80036
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9784452 0.0568455 34.80 <2e-16 ***
## Determinacion 0.0139704 0.0009444 14.79 <2e-16 ***
## Determinacion:Heridas3-4 -0.0064198 0.0006182 -10.38 <2e-16 ***
## Determinacion:Heridas5+ -0.0107028 0.0006843 -15.64 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2086805)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 103.41 on 807 degrees of freedom
## AIC: 4554.2
##
## Number of Fisher Scoring iterations: 5De entre estos cuatro modelos, el que mejor ajuste produce es el modelo completo con interacciones, con un AIC de 4544.8, el más bajo de entre los cuatro. Al igual que ocurría con los modelos con la función link inversa, estos modelos tienen los mismos problemas de aplicabilidad. Por temas de espacio, solo se incluyen las gráficas del primer modelo, las del resto son muy similares así que no merece la pena incluirlas.
Ahora vamos a probar a introducir términos polinómicos de la variable Determinacion con tal de ver si mejora tanto el ajuste como la aplicabilidad. Vamos a automatizar la búsqueda del mejor modelo con términos polinómicos, de tal forma que nos quedaremos con aquel que tenga menor AIC.
AICs <- c()
for(i in 1:6){
model <- glm(T_batalla ~ poly(Determinacion,i), data = datos, family = Gamma(link = "log"))
AICs[i] <- model$aic
}
AICs## [1] 4896.281 4898.200 4900.180 4902.114 4898.507 4898.648El modelo con el solo el término polinómico de grado 1 es el que menor AIC tiene. Vamos a continuar la modelización teniendo eso en cuenta. Vamos a probar también con la transformación logaritmica de Determinación y ver si mejora el modelo.
summary(mod.log <- glm(T_batalla ~ log(Determinacion), data = datos, family = Gamma(link = "log")))##
## Call:
## glm(formula = T_batalla ~ log(Determinacion), family = Gamma(link = "log"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.67185 -0.31666 -0.14300 0.06312 2.90685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.46667 0.25450 1.834 0.0671 .
## log(Determinacion) 0.50295 0.06271 8.021 3.68e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.3835453)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 159.95 on 809 degrees of freedom
## AIC: 4913.3
##
## Number of Fisher Scoring iterations: 6El logaritmo no mejora el ajuste con respecto al mismo modelo sin la transformación, así que finalmente, vamos a ver el error de prediccion cometido por el modelo completo con interacciones mediante validación cruzada.
set.seed(1)
cv.glm(datos, log.interac, K = 10)$delta## [1] 53.66316 53.60299Haber empleado la función link log no ha supuesto ninguna mejora ni ningún empeoramiento sobre el error medio del modelo con la función link inversa. Determinamos que proporción de la DEVIANCE explica el modelo.
1-(log.interac$deviance/log.interac$null.deviance)## [1] 0.4528238Por lo tanto, para este modelo en el cual se ha empleado el enlace inverse, se consigue un AIC=4544 y un porcentaje de la DEVIANCE explicada del 45.28%.
Ajuste con función link identidad
Para este ajuste vamos a seguir un proceso idéntico al seguido con la función link inversa y log. Comenzamos ajustando 4 modelos: un modelo reducido, un modelo completo, un modelo con las interacciones y un modelo reducido con interacción.
ident.reducido <- glm(T_batalla ~ Determinacion, data = datos, family = Gamma(link = "identity"))
summary(ident.reducido)##
## Call:
## glm(formula = T_batalla ~ Determinacion, family = Gamma(link = "identity"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.65718 -0.31483 -0.13358 0.07065 2.91447
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.01707 0.78389 6.400 2.63e-10 ***
## Determinacion 0.12235 0.01393 8.784 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.3767992)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 157.42 on 809 degrees of freedom
## AIC: 4900
##
## Number of Fisher Scoring iterations: 5par(mfrow = c(2,2), mar = c(2.1,4,2.1,4))
plot(ident.reducido)
ident.completo <- glm(T_batalla ~ Determinacion + Heridas,
data = datos, family = Gamma(link = "identity"))
summary(ident.completo)##
## Call:
## glm(formula = T_batalla ~ Determinacion + Heridas, family = Gamma(link = "identity"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.47724 -0.26107 -0.07695 0.11362 2.75346
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.86488 0.85842 14.987 < 2e-16 ***
## Determinacion 0.07918 0.01017 7.788 2.1e-14 ***
## Heridas3-4 -6.03750 0.65036 -9.283 < 2e-16 ***
## Heridas5+ -8.78608 0.64652 -13.590 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2221883)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 106.69 on 807 degrees of freedom
## AIC: 4580
##
## Number of Fisher Scoring iterations: 6ident.interac <- glm(T_batalla ~ Determinacion*Heridas,
data = datos, family = Gamma(link = "identity"))
summary(ident.interac)##
## Call:
## glm(formula = T_batalla ~ Determinacion * Heridas, family = Gamma(link = "identity"),
## data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4867 -0.2503 -0.0706 0.1135 2.7273
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.74581 1.66081 2.858 0.00438 **
## Determinacion 0.21059 0.02810 7.494 1.76e-13 ***
## Heridas3-4 2.79533 1.87119 1.494 0.13560
## Heridas5+ 0.29284 1.84945 0.158 0.87423
## Determinacion:Heridas3-4 -0.14463 0.03187 -4.539 6.52e-06 ***
## Determinacion:Heridas5+ -0.14951 0.03149 -4.747 2.44e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2036734)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 101.86 on 805 degrees of freedom
## AIC: 4545.7
##
## Number of Fisher Scoring iterations: 6ident.red.int <- glm(T_batalla ~ Determinacion + Determinacion:Heridas,
data = datos, family = Gamma(link = "identity"))
summary(ident.red.int)##
## Call:
## glm(formula = T_batalla ~ Determinacion + Determinacion:Heridas,
## family = Gamma(link = "identity"), data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.50631 -0.25729 -0.07007 0.10687 2.79981
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.12186 0.56554 10.825 <2e-16 ***
## Determinacion 0.18898 0.01278 14.787 <2e-16 ***
## Determinacion:Heridas3-4 -0.09893 0.01027 -9.631 <2e-16 ***
## Determinacion:Heridas5+ -0.14554 0.01015 -14.341 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2082567)
##
## Null deviance: 185.96 on 810 degrees of freedom
## Residual deviance: 102.93 on 807 degrees of freedom
## AIC: 4550.3
##
## Number of Fisher Scoring iterations: 5De entre estos cuatro modelos, el que mejor ajuste produce es el modelo completo con interacciones, con un AIC de 4545.7, el más bajo de entre los cuatro. Al igual que ocurría con los modelos con la función link inversa, estos modelos tienen los mismos problemas de aplicabilidad. Las gráficas de todos los modelos son muy similares, así que solamente se ha incluído las del primer modelo.
Ahora vamos a probar a introducir términos polinómicos de la variable Determinacion con tal de ver si mejora tanto el ajuste como la aplicabilidad. Vamos a automatizar la búsqueda del mejor modelo con términos polinómicos, de tal forma que nos quedaremos con aquel que tenga menor AIC.
AICs <- c()
for(i in 1:6){
model <- glm(T_batalla ~ poly(Determinacion,i), data = datos,
family = Gamma(link = "identity"))
AICs[i] <- model$aic
}
AICs## [1] 4899.976 4898.185 4900.168 4902.024 4900.523 4899.736Atendiendo a los valores del AIC el mejor modelo sería el del polinomio de grado 2, sin embargo, la diferencia con el modelo de grado 1 es insignificante, así que por simplicidad, nos quedamos con el modelo del polinomio de grado 1.
La transformación logaritmica de Determinación con la función link identity no produce un modelo porque no es capaz de encontrar un conjunto válido de coeficientes. Por lo que no vamos a contar con dicha transformación. Finalmente, vamos a ver el error de predicción cometido por el modelo completo con interacciones mediante validación cruzada.
set.seed(1)
cv.glm(datos, ident.interac, K = 10)$delta## [1] 53.50928 53.46274Haber empleado la función link log no ha supuesto ninguna mejora ni ningún empeoramiento sobre el error medio de los modelo con las funciones link inversa y log. Determinamos que proporción de la DEVIANCE explica el modelo.
1-(ident.interac$deviance/ident.interac$null.deviance)## [1] 0.4522536Por lo tanto, para este modelo en el cual se ha empleado el enlace inverse, se consigue un AIC=4544 y un porcentaje de la DEVIANCE explicada del 45.22%.
Comparación de modelos: Validación del ajuste y Capacidad predictiva
A continuación se muestra una tabla comparativa a nivel ajuste y predicción de los diferente modelos.
| Link | Modelo | %D explicada | AIC | RMSE |
|---|---|---|---|---|
| Inverso | Completo | 44.67 | 4549.9 | 53.6 |
| Logaritmo | Interacción | 45.28 | 4544.8 | 53.6 |
| Identidad | Interacción | 45.22 | 4545.7 | 53.5 |
Por lo tanto, podemos determinar que no hay ningún modelo mucho mejor que otro en términos de ajuste y predicción, siendo el modelo con la función link logarítmica la que presenta, aunque por poco, un valor de AIC y DEVIANCE menor. Además, es el que mejor interpretabilidad tiene.
Finalmente, podemos ver graficamente los valores ajustados y observados del modelo seleccionado.
par(mfrow=c(1,1))
plot(log.interac$fitted.values, datos$T_batalla, main='Mejor modelo', xlab='Valores Ajustados', ylab='Valores Observados')
abline(0,1)
Modelo final
Por ello, elegimos este modelo como modelo final. \[ Y_{ij} \sim Ga(\nu_{ij}, \lambda_{ij}) \] \[ log(\mu_{ij}) = 2.06 + (0.012 -0.007d_{(3-4)}-0.005d_{(5+)})x_{i} + 0.03d_{(3-4)} -0.33 d_{(5+)} \]
con \(X_i\) como variable Determinacion y \(d_j\) como variable indicadora de Heridas.
Así pues: \[ \widehat{E(Y_{ij})} = \mu_{ij} = exp(2.06 + (0.012 -0.007d_{(3-4)}-0.005d_{(5+)})x_{i} + 0.03d_{(3-4)} -0.33 d_{(5+)}) \] Por lo tanto, el modelo se interpreta como cambios porcentuales en \(E(Y_{ij})\) por cada incremento de una unidad en \(X_i\) (% = \(100 \cdot (e^{\beta_i})\)).
Análisis del BMI mediante regresión gamma
Introducción
El siguiente análisis consiste en la modelización mediante regresión gamma de los datos de la American Time Use Survey (ATUS) Eating & Health Module de 2014. Concretamente, se modelizará el BMI en función del resto de variables.
Se puede encontrar más información sobre los datos en su pagína web, aunque el dataset empleado en este caso es una conversión a CSV que se puede encontrar en Kaggle.
En primer lugar, cargamos los datos (eliminando aquellas observaciones para las que el BMI reportado sea 0, ya que es imposible):
bmi_data <- read_csv("data/ehresp_2014.csv", show_col_types = FALSE) %>%
filter(erbmi > 0)Análisis descriptivo
bmi <- bmi_data
str(bmi, give.attr=FALSE)## spec_tbl_df [10,637 × 37] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ tucaseid : num [1:10637] 2.01e+13 2.01e+13 2.01e+13 2.01e+13 2.01e+13 ...
## $ tulineno : num [1:10637] 1 1 1 1 1 1 1 1 1 1 ...
## $ eeincome1 : num [1:10637] -2 1 2 2 1 1 1 1 1 1 ...
## $ erbmi : num [1:10637] 33.2 22.7 49.4 31 30.7 ...
## $ erhhch : num [1:10637] 1 3 3 3 3 1 3 3 3 3 ...
## $ erincome : num [1:10637] -1 1 5 5 1 1 1 1 1 1 ...
## $ erspemch : num [1:10637] -1 -1 -1 -1 1 5 -1 -1 5 -1 ...
## $ ertpreat : num [1:10637] 30 45 60 65 20 30 30 117 80 35 ...
## $ ertseat : num [1:10637] 2 14 0 0 10 5 5 10 0 20 ...
## $ ethgt : num [1:10637] 0 0 0 0 0 0 0 0 0 0 ...
## $ etwgt : num [1:10637] 0 0 0 0 0 0 0 0 0 0 ...
## $ eudietsoda : num [1:10637] -1 -1 -1 -1 1 -1 -1 -1 2 1 ...
## $ eudrink : num [1:10637] 2 2 1 1 1 1 2 2 1 1 ...
## $ eueat : num [1:10637] 1 1 2 2 1 1 1 1 2 1 ...
## $ euexercise : num [1:10637] 2 2 2 1 1 2 1 1 2 2 ...
## $ euexfreq : num [1:10637] -1 -1 -1 5 2 -1 3 6 -1 -1 ...
## $ eufastfd : num [1:10637] 2 1 2 2 1 1 1 2 1 1 ...
## $ eufastfdfrq: num [1:10637] -1 1 -1 -1 3 3 1 -1 2 5 ...
## $ euffyday : num [1:10637] -1 2 -1 -1 1 2 2 -1 1 2 ...
## $ eufdsit : num [1:10637] 1 1 1 1 1 1 1 1 1 1 ...
## $ eufinlwgt : num [1:10637] 5202086 29400000 26000000 17500000 3661280 ...
## $ eusnap : num [1:10637] 1 2 2 1 2 2 2 2 2 2 ...
## $ eugenhth : num [1:10637] 1 2 5 4 3 2 2 3 1 4 ...
## $ eugroshp : num [1:10637] 1 3 2 1 2 3 1 1 1 1 ...
## $ euhgt : num [1:10637] 60 63 62 69 71 65 63 70 65 70 ...
## $ euinclvl : num [1:10637] 5 5 5 5 5 5 5 5 5 5 ...
## $ euincome2 : num [1:10637] -2 -1 2 2 -1 -1 -1 -1 -1 -1 ...
## $ eumeat : num [1:10637] 1 1 -1 1 -1 1 1 1 1 -1 ...
## $ eumilk : num [1:10637] 2 2 -1 2 -1 2 2 2 2 -1 ...
## $ euprpmel : num [1:10637] 1 1 2 1 2 3 1 1 1 2 ...
## $ eusoda : num [1:10637] -1 -1 2 2 1 2 -1 -1 1 1 ...
## $ eustores : num [1:10637] 2 1 -1 1 -1 2 1 1 3 1 ...
## $ eustreason : num [1:10637] 1 2 -1 1 -1 5 3 4 1 1 ...
## $ eutherm : num [1:10637] 2 2 -1 2 -1 2 2 2 2 -1 ...
## $ euwgt : num [1:10637] 170 128 270 210 220 200 155 180 170 282 ...
## $ euwic : num [1:10637] 1 2 2 1 2 2 -1 -1 -1 -1 ...
## $ exincome1 : num [1:10637] 2 0 12 0 0 0 0 0 0 0 ...Tenemos un banco de datos con 10637 observaciones y 37 variables. Las dos primeras variables son variables de identificación, y la última es de etiquetado, por lo que las eliminamos por simplicidad.
bmi <- bmi[, -c(1, 2, 37)]Respecto a las variables restantes, podemos destacar que:
- “erbmi” es una variable cuantitativa continua, y constituye nuestra variable respuesta.
- “ertpreat”, “ertseat”, “euexfreq”, “eufastfdfrq”, “eufinlwgt”, “euhgt”, “euwgt” son variables cuantitativas.
- El resto de variables son cualitativas, por lo que han de ser convertidas a factor.
En primer lugar, observamos las correlaciones entre las variables.
corrplot(cor(bmi), type = "upper")
Se aprecian correlaciones entre las variables que aportan información similar. En el caso del BMI una correlación importante con el peso. Respecto a las variables cuantitativas:
no.factor <- c("erbmi", "ertpreat", "ertseat", "euexfreq", "eufastfdfrq",
"eufinlwgt", "euhgt", "euwgt")
summary(bmi[,no.factor])## erbmi ertpreat ertseat euexfreq
## Min. :13.00 Min. : 0.00 Min. : -2.0 Min. :-2.00
## 1st Qu.:23.60 1st Qu.: 30.00 1st Qu.: 0.0 1st Qu.:-1.00
## Median :26.60 Median : 60.00 Median : 4.0 Median : 2.00
## Mean :27.77 Mean : 65.85 Mean : 16.9 Mean : 2.27
## 3rd Qu.:30.70 3rd Qu.: 90.00 3rd Qu.: 15.0 3rd Qu.: 4.00
## Max. :73.60 Max. :508.00 Max. :990.0 Max. :38.00
## eufastfdfrq eufinlwgt euhgt euwgt
## Min. :-2.000 Min. : 756844 Min. :56.00 Min. : 98.0
## 1st Qu.:-1.000 1st Qu.: 3511427 1st Qu.:64.00 1st Qu.:145.0
## Median : 1.000 Median : 6023405 Median :66.00 Median :170.0
## Mean : 1.158 Mean : 8223308 Mean :66.69 Mean :176.3
## 3rd Qu.: 2.000 3rd Qu.: 10300000 3rd Qu.:70.00 3rd Qu.:200.0
## Max. :21.000 Max. :103000000 Max. :77.00 Max. :340.0par(mfrow=c(2,4))
boxplot(bmi$erbmi, main="erbmi")
boxplot(bmi$ertpreat, main="ertpreat")
boxplot(bmi$ertseat, main="ertseat")
boxplot(bmi$euexfreq, main="euexfreq")
boxplot(bmi$eufastfdfrq, main="eufastfdfrq")
boxplot(bmi$eufinlwgt, main="eufinlwgt")
boxplot(bmi$euhgt, main="euhgt")
boxplot(bmi$euwgt, main="euwgt")
Vemos que todas las variables muestran valores atípicos. En el caso del BMI es esperable por su distribución. Otras variables, como ertseat tienen un número mucho mayor, por lo que debe tenerse en cuenta al modelizar y diagnosticar el modelo.
Una vez completada la descripción de las variables, convertimos a factor las variables categóricas:
factor.names <- c("erhhch", "erspemch", "ethgt", "etwgt", "eudietsoda", "eudrink",
"eueat", "euexercise", "eufastfd", "euffyday", "eugroshp",
"euinclvl", "eusnap", "eumeat", "eumilk", "euprpmel", "eusoda",
"eustores", "eustreason", "eutherm", "euwic",
"eeincome1", "erincome", "eufdsit", "eugenhth", "euincome2")
bmi[,factor.names] <- lapply(bmi[,factor.names], factor)
summary(bmi[,factor.names])## erhhch erspemch ethgt etwgt eudietsoda eudrink eueat
## 1: 500 -1:5249 0:10555 0:10535 -3: 1 -3: 1 -2: 55
## 2: 213 1 : 221 1: 40 1: 53 -2: 4 -2: 7 1 :5822
## 3:9924 2 : 90 2: 42 2: 49 -1:7767 1 :7164 2 :4760
## 3 : 224 1 :1121 2 :3465
## 4 : 161 2 :1677
## 5 :4692 3 : 67
##
## euexercise eufastfd euffyday eugroshp euinclvl eusnap eumeat
## -2: 6 -2: 23 -2: 2 -3: 1 5:8760 -2: 35 -2: 9
## 1 :6715 1 :6196 -1:4441 -2: 2 6:1877 1 :1091 -1:2973
## 2 :3916 2 :4418 1 :2261 1 :6498 2 :9511 1 :6784
## 2 :3933 2 :2841 2 : 871
## 3 :1295
##
##
## eumilk euprpmel eusoda eustores eustreason eutherm euwic
## -3: 1 -2: 3 -2: 3 -2: 50 2 :2895 -2: 3 -2: 22
## -2: 1 1 :6603 -1:3473 -1:2841 -1 :2891 -1:3853 -1:5130
## -1:2973 2 :2970 1 :2870 1 :5254 1 :2486 1 : 807 1 : 356
## 1 : 150 3 :1061 2 :4291 2 :1927 3 :1030 2 :5974 2 :5129
## 2 :7512 3 : 337 4 : 668
## 4 : 35 6 : 435
## 5 : 193 (Other): 232
## eeincome1 erincome eufdsit eugenhth euincome2
## -3: 102 -1: 213 -3: 1 -3: 2 -3: 197
## -2: 144 1 :6696 -2: 10 -2: 28 -2: 538
## 1 :6696 2 : 504 1 :9987 1 :1943 -1:6577
## 2 :3258 3 : 932 2 : 514 2 :3595 1 :1066
## 3 : 437 4 : 33 3 : 125 3 :3288 2 :1917
## 5 :2259 4 :1303 3 : 342
## 5 : 478Búsqueda del modelo óptimo
Como hemos visto, el BMI se distribuye aproximadamente como una gamma, por lo que aplicamos este tipo de regresión. EM primer lugar, obtenemos modelos completos con todos los links disponibles, para luego hacer una selección de modelos.
full.id <- glm(erbmi ~ .,
data = bmi,
family=Gamma(link="identity"))
par(mfrow=c(2,2)); plot(full.id); par(mfrow=c(1,1))
plot(residuals(full.id, type = "deviance"))
full.log <- glm(erbmi ~ .,
data = bmi,
family=Gamma(link="log"))
par(mfrow=c(2,2)); plot(full.log); par(mfrow=c(1,1))
plot(residuals(full.log, type = "deviance"))
full.inv <- glm(erbmi ~ ., # forumla
data = bmi, # dataset
family=Gamma(link="inverse"))
par(mfrow=c(2,2)); plot(full.inv); par(mfrow=c(1,1))
plot(residuals(full.inv, type = "deviance"))
Como vemos, los modelos completos cumplen condiciones como valores de los residuos DEVIANCE entre -2 y 2 en los tres link. Respecto al resto de diagnósticos, parece que el link identidad tiene mejores residuos.
Ahora seleccionamos modelos, para lo que hacemos una búsqueda mediante step. Usaremos el mejor modelo estimado por la función.
step.id.out <- step(full.id, direction = "backward")
step.log.out <- step(full.log, direction = "backward")
step.inv.out <- step(full.inv, direction = "backward")Por tanto, los modelos que ajustaremos serán:
step.id.out$call## glm(formula = erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth +
## eugroshp + euhgt + eusoda + eustreason + euwgt, family = Gamma(link = "identity"),
## data = bmi)step.inv.out$call## glm(formula = erbmi ~ eeincome1 + erhhch + erincome + erspemch +
## ethgt + etwgt + eueat + euexercise + euexfreq + eufastfd +
## eufinlwgt + eusnap + eugenhth + euhgt + euincome2 + eumeat +
## euprpmel + eusoda + euwgt + euwic, family = Gamma(link = "inverse"),
## data = bmi)step.log.out$call## glm(formula = erbmi ~ erspemch + ethgt + etwgt + eueat + euexfreq +
## eufinlwgt + eugenhth + euhgt + euincome2 + eusoda + euwgt +
## euwic, family = Gamma(link = "log"), data = bmi)step.id <- glm(erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth +
eugroshp + euhgt + eusoda + eustreason + euwgt,
family = Gamma(link = "identity"),
data = bmi)
par(mfrow=c(2,2)); plot(step.id); par(mfrow=c(1,1))
plot(residuals(step.id, type = "deviance"))
step.log <- glm(formula = erbmi ~ erspemch + ethgt + etwgt + eueat + euexfreq +
eufinlwgt + eugenhth + euhgt + euincome2 + eusoda + euwgt +
euwic, family = Gamma(link = "log"), data = bmi)
par(mfrow=c(2,2)); plot(step.log); par(mfrow=c(1,1))
plot(residuals(step.log, type = "deviance"))
step.inv <- glm(formula = erbmi ~ eeincome1 + erhhch + erincome + erspemch +
ethgt + etwgt + eueat + euexercise + euexfreq + eufastfd +
eufinlwgt + eusnap + eugenhth + euhgt + euincome2 + eumeat +
euprpmel + eusoda + euwgt + euwic, family = Gamma(link = "inverse"),
data = bmi)
par(mfrow=c(2,2)); plot(step.inv); par(mfrow=c(1,1))
plot(residuals(step.inv, type = "deviance"))
En términos de AIC, el mejor modelo es el que emplea el enlace identidad, con AIC = 21647.03, contra el modelo con enlace logaritmo (AIC = 28302.22) y el de enlace inverso (AIC = 40643.39).
Si lo comparamos en cuanto a varianza explicada:
100 * (1 - (step.id$deviance/step.id$null.deviance))## [1] 98.65191100 * (1 - (step.log$deviance/step.log$null.deviance))## [1] 97.47904100 * (1 - (step.inv$deviance/step.inv$null.deviance))## [1] 91.98856De nuevo, el mejor modelo es el que emplea en enlace identidad, que explica un 98.65% de la DEVIANCE. Si observamos sus coeficientes:
summary(step.id)##
## Call:
## glm(formula = erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth +
## eugroshp + euhgt + eusoda + eustreason + euwgt, family = Gamma(link = "identity"),
## data = bmi)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.131176 -0.010936 0.001528 0.009796 0.229545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.4083315 0.8729159 58.893 < 2e-16 ***
## erincome1 0.0745781 0.0443039 1.683 0.092341 .
## erincome2 0.1312195 0.0521158 2.518 0.011822 *
## erincome3 0.0839181 0.0484358 1.733 0.083203 .
## erincome4 0.1090085 0.1143751 0.953 0.340572
## erincome5 0.1567179 0.0461082 3.399 0.000679 ***
## ethgt1 0.8722284 0.0924147 9.438 < 2e-16 ***
## ethgt2 0.9540602 0.1216994 7.839 4.96e-15 ***
## etwgt1 -1.4951243 0.1708827 -8.749 < 2e-16 ***
## etwgt2 -0.9310082 0.0646855 -14.393 < 2e-16 ***
## eusnap1 0.0017920 0.1044026 0.017 0.986306
## eusnap2 -0.0698268 0.1021802 -0.683 0.494389
## eugenhth-2 0.2271269 0.4625373 0.491 0.623404
## eugenhth1 0.3133602 0.4455386 0.703 0.481866
## eugenhth2 0.3217543 0.4454499 0.722 0.470118
## eugenhth3 0.3542196 0.4454281 0.795 0.426495
## eugenhth4 0.3976043 0.4456215 0.892 0.372281
## eugenhth5 0.4523042 0.4461883 1.014 0.310746
## eugroshp-2 0.5676538 1.0056308 0.564 0.572443
## eugroshp1 0.7083560 0.8848995 0.800 0.423443
## eugroshp2 0.7115029 0.8897742 0.800 0.423935
## eugroshp3 0.6458298 0.8848067 0.730 0.465461
## euhgt -0.7684382 0.0019201 -400.199 < 2e-16 ***
## eusoda-1 -0.1218794 0.3278954 -0.372 0.710122
## eusoda1 -0.1318953 0.3279686 -0.402 0.687576
## eusoda2 -0.1569894 0.3278840 -0.479 0.632094
## eustreason-2 -0.7734256 0.5980125 -1.293 0.195926
## eustreason-1 -0.8762340 0.5992355 -1.462 0.143701
## eustreason1 -0.8137849 0.5921778 -1.374 0.169401
## eustreason2 -0.8602364 0.5921865 -1.453 0.146352
## eustreason3 -0.8598634 0.5924067 -1.451 0.146677
## eustreason4 -0.7647428 0.5926035 -1.290 0.196912
## eustreason5 -0.7910744 0.5942001 -1.331 0.183110
## eustreason6 -0.8720373 0.5928597 -1.471 0.141348
## euwgt 0.1560453 0.0002035 766.798 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.0006107502)
##
## Null deviance: 476.2832 on 10636 degrees of freedom
## Residual deviance: 6.4207 on 10602 degrees of freedom
## AIC: 21647
##
## Number of Fisher Scoring iterations: 5Vemos que solo los coeficientes correspondientes a las variables erincome, ethgt, euhgt y euwgt, que se refieren a ingresos, peso y altura, son significativos.
Al emplear un enlace identidad, podemos interpretar el modelo como que cada incremento de una unidad en \(X_j\) hace crecer el valor esperado del BMI, \(E(BMI)\) en \(\beta_j\), asumiendo que el resto de variables se mantienen constantes.
Comparación con un modelo lineal
A continuación, ajustaremos el mismo modelo que hemos estimado como óptimo para ver si hay diferencias.
gaussian <- glm(erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth +
eugroshp + euhgt + eusoda + eustreason + euwgt,
family = ("gaussian"),
data = bmi)
par(mfrow=c(2,2)); plot(gaussian); par(mfrow=c(1,1))
Los residuos no se ajustan a la normalidad, además de ser heterocedásticos. Por tanto, vemos que la elección de un modelo gamma es más adecuado a la hora de describir los datos.
Validación cruzada
Pese a que el objetivo del modelo no es la predicción, realizamos una verificación mediante validación cruzada para comprobar como se ajustan los valores.
set.seed(1)
data_partition <- createDataPartition(bmi$erbmi,
times = 1, p = 0.8,
list = FALSE)
train <- bmi[data_partition,]
test <- bmi[-data_partition,]
cv.model <- glm(erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth +
eugroshp + euhgt + eusoda + eustreason + euwgt,
family = Gamma(link = "identity"),
data = test)
prediction <- predict(cv.model, newdata=test, type="response")
plot(test$erbmi, prediction)
Vemos que la relación es prácticamente lineal, aunque los valores extremos tienen más variabilidad.